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_Ie-  describes  hex^>  the  use  of  continuation  or  homotopy  procedures  in  an 
efficient  algorithm  for  the  calculation  of  raypaths  froa  points  within  a 
region  at  depth  to  receivers  on  the  surface.;.  Our  work  follows  closely  that 
of  H.  B.  Keller  and  his  students  D.  J.  Perozzi  [1983]  and  J.  A.  Fawcett 
[1983].  Polynoaials  define  the  interfaces  in  a  two-dimensional  piecewise 
constant  velocity  aediua.  Starting  with  horizontal  layering  only,  the 
interfaces  are  gradually  deformed  until  the  desired  earth  smdel  is  achieved. 
At  each  deformation  step  the  ray  equations  determined  by  Fermat's  Principle 
are  solved  using  Newton's  method  and  the  ray  from  the  first  source  position 

to  a  receiver  vertically  above  it  is  found. 

5  <j  W'O  r 

Next  wa  employ  ^source  continuation,  moving  the  source  on  a  grid  within 

/«. 

the  region  of  interest.  At  each  source  position  we  find^the  rays  to  all 
receivers  using  continuation  in  receiver  location. 

Knowing  the  raypaths,  it  is  straightforward  to  construct  a  table  of 
traveltimes.  These  traveltimes  alone  serve  to  position  reflectors  in  the 
subsurface. 

Ve  giver  two  examples  of  migration  and  indicate  possible  applications  to 


an 


1 .  INTRO  MICTION 


Structural  Inversion  or  aigration  is  possible  if  traveltimes  between 
subsurface  points  and  surface  receivers  are  known.  Indeed,  travel  tines 
alone  serve  to  position  reflectors,  whereas  the  strength  of  a  reflection 
contains  information  concerning  the  physical  properties  of  the  aediua.  The 
calculation  of  traveltiae  is  straightforward  if  the  sonnd  path  between 
source  and  receiver  is  known.  To  perform  a  migration  many  such  sound  paths 
must  be  found  and  in  the  past  this  has  proved  too  expensive  to  make  the 
method  practical.  Ve  present  a  high  speed  ray  tracing  routine  at  the  heart 
of  a  very  accurate  aigration  procedure. 


Our  ray  tracing  routine  is  based  upon  continuation  or  hoaotopy 
procedures.  The  method  was  developed  by  H.  B.  Keller  and  two  of  his 
students.  D.  J.  Peroxzi  [1983]  and  J.  A.  Fawcett  [1983].  It  uses  a  known 
ray  solution  to  provide  a  guess  at  the  solution  for  a  nearby  raypath.  With 
aigration  in  mind,  we  have  extended  their  ideas  and  created  an  algorithm 
suitable  for  finding  large  numbers  of  raypaths  and  traveltiaes. 


Some  prior  knowledge  of  the  subsurface  is  assumed  and  the  user  provides 
a  model  of  known  interfaces  and  interval  velocities.  A  region  of  interest 
is  chosen  and  we  consider  every  point  within  the  region  as  a  possible  point 
source.  Ray  tracing  begins  between  these  points  at  depth  and  receivers  on 
the  surface. 


The  first  ray  is  found  in  a  stratified  medium,  where  the  depth  to  an 
interface  is  some  average  depth  in  the  true  medium.  The  interfaces  are  then 


gradually  deforaed  until  they  become  thoae  in  the  trne  aodel.  At  each 
do formation  step  we  solve  a  systea  of  equations  deterained  by  Fermat's 
Principle .  Next  we  use  a  continuation  procedure  for  the  source  location  to 
aove  the  source  within  the  region  of  interest.  At  each  source  position  we 
find  raypaths  to  all  receivers  using  continuation  in  receiver  location. 

For  each  raypath  a  traveltiae  is  calculated  and  output.  Thus,  for  each 
source  position,  there  is  a  list  of  traveltiaes  to  receivers.  A  aigrated 
depth  section  is  produced  when  coaaon  depth  point  data  is  sunaed  over  curves 
deterained  by  these  traveltiaes. 

Ve  give  two  exaaples  of  aigration  of  synthetic  data.  They  show  that 
the  aethod  positions  reflectors  accurately,  though  we  sake  no  claias 
regarding  the  aaplitude  of  the  output.  Our  ray  tracing  procedure  is  also  a 
tool  for  forward  aodeling.  Ve  indicate  several  possible  applications. 


% 


2.  FORMULATION  OF  THE  NUMERICAL  METHOD 

The  aediiia  is  coaprised  of  constant  velocity  layers  separated  by 
arbitrary  interfaces.  Each  interface  is  characterized  by  a  polynoaial,  f. 
and  a  constant,  c,  in  an  equation  of  the  fora 


z  -  ck  +  X  fk<*>  •  (1) 

The  significance  of  the  parameter  X  will  becoae  apparent  later.  Ve 
seek  the  ray  joining  a  source  at  depth  and  a  receiver  on  the  snrfaee.  The 
ray  is  a  series  of  straight  line  segaents  coapletely  defined  by  ita  points 
of  intersection  with  each  interface,  given  soae  sonrce  and  receiver. 

Let  the  nnaber  of  interfaces  between  source  and  receiver  be  N.  Let  the 
coordinates  of  the  point  of  intersection  of  the  ray  with  the  ktjj  interface 
be  *  (xfc,  Zfc)  ■  (Xfc,  Cfc  +  Xffc(xfc))  ,  and  the  velocity  of  the  aediua  on 
the  ktj,  segaent  be  Vr.  The  source  position  is  Xf  *  Ixq,  z (j)»  the 

receiver  position  is  Xr  =  (xN+1,  *N+1>*  <See  Fil«r«  1). 


Define 


Using  (2) 


V  V 

k+1  *  k 

-*-±-  (Ax.  +  XAx.f  (x  ))  - 
k  *  n  x  "k+1 


(Axk+1  +  *Azk+1  fk(xk>)  “  0 


Thus,  there  sre  N  equations  in  the  N  unknowns,  xk»  and  in  systea 
wri  te 


f(XU).k)  -  0  . 


where. 


I  "  <#•*. 


Not  surprisingly  the  suae  systea  of  equations  is  obtained  if 
that  Snell's  Law  be  satisfied  at  each  interface. 


i  .  <«) 

fora  we  aay 

(7) 

(8) 

we  require 
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FINDING  THE  FIRST  RAT  -  CONTINUATION  IN  INTERFACES 


Setting  X  *  0  in  (1)  nil  the  interfaces  becoae  horizontal)  the  aediua 
ia  stratified.  The  raypath  between  the  source  and  a  receiver  vertically 
above  it  on  the  surface  is  fonnd  (it  ia  a  straight  vertical  line).  The 
paraaeter  X  is  then  increaented  by  soae  AX  and  we  seek  the  solution  of  (7) 
using  Newton's  aethod.  Assuaing  that  we  know  a  solution  X(X),  for  soae 
value  of  X,  (for  ezaaple.  we  know  it  for  X  -  0),  then  as  an  initial  guess  in 
Newton's  aethod  for  the  solution  at  X  +  AX  we  use 


X(0)(X  +  A  X)  *  X  (X)  +  AX  X  (X)  • 


(9) 


In  this  equation  X(X)  *  dX/dX  .  We  differentiate  (7)  to  obtain  an  equation 
froa  which  1(X)  can  be  deterained: 


*!  »f  .  9f 

■ar  "  '3T  *a)  +  ax  “  0  ’  (10) 


Now  using  Newton's  aethod,  an  iaproved  value  for  the  solution  at  X  +  AX  is 


X 


(1) 


(11) 


In  this  equation,  J  ■  }|/8X  is  the  Jacobian  of  the  systen.  Froa  (6)  we  see 
that  k  depends  only  on  ife-i,  xv  and  Xfc+i>  the  Jacobian  is  tridiagonal.  To 


find  its  inverse  ve  nee  a  banded  system  solver  from  the  I.N.S.L.  library. 
If  Newton's  method  fails  to  converge  to  a  solution  within,  say,  four 
iterations  we  set  AX  ■  AX/2  and  return  to  (9).  Ve  proceed  in  this  manner 
until  X  ■  1,  when  the  first  ray  in  the  true  model  is  found.  As  Fawcett 
points  out.  in  almost  all  cases  it  is  possible  to  begin  by  setting  X  *  0, 
and  AX  *  1  and  in  one  step  obtain  an  initial  guess  for  which  Newton's  method 
converges.  The  algorithm  is  indeed  extremely  efficient. 


4.  CONTINUATION  IN  BECBIVER  LOCATION 


Given  the  ray  eolation  for  a  receiver  at  we  eaploy  coat  innation 
procedures  to  find  the  ray  solution  for  an  adjacent  receiver  at  1|^^  .  For 
the  reeeiver  position  ve  write 


£r<x) 


+  (1-X) 


(12) 


As  X  varies  between  0  and  1.  Xg^  is  continued  into  Xg^«  (See  Fig.  2). 
Onr  systea  is  now 


f(  X(X).  X^X))  =  0  • 


(13) 


Again  we  use  (9)  along  with  Newton's  aethod.  This  tiae  i(X)  is  deterained 
as  a  solution  of  the  derivative  of  (13): 


and 


Since  only  ^  depends  on  Xjj,  the  Nx2  matrix  df/dlg  has  all  elements  equal  to 
zero  except  for  dp^/dxjj+i  and  ^  in  the  last  row.  Here  we  consider 

the  surface  of  the  earth  to  be  flat  (z  =  0) .  Thus.  dp^/3zj,j+1  =  0  also.  The 
non-zero  element  in  (15)  is  the  receiver  spacing.  Ve  could  have  included 
surface  topography  by  writing 

1  =  CN+1  +  XWX)  ‘ 


In  finding  the  first  ray  we  always  place  the  receiver  vertically  above  the 
source.  This  allows  us  to  write  down  the  solution  in  a  stratified  medium. 
In  general  the  receiver  will  be  at  some  other  position  and  we  simply  use 
continuation  in  receiver  location  to  proceed  to  that  point. 

Continuation  in  receiver  location  was  the  work  of  Perozzi  and 
continuation  in  interfaces  was  Pawcett’s  idea.  Our  aim  here  is  to  create  a 
table  of  traveltimes  between  points  at  depth  and  receivers  on  the  surface. 
Ve  could  begin  at  each  source  position  with  a  stratified  earth  and  proceed 
as  above.  More  efficiently,  we  could  apply  the  same  ideas  of  continuation 
to  the  source. 
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CONTINUATION  IN  SflBKCB  LOCATION 


If  the  ray  eolation  for  a  source  at  Xg'1'  is  known,  we  can  employ 
continuation  procedures  to  find  the  solution  for  a  nearby  source  location 
Writing 


XgU)  =  X  X^j)  +  (1  -  X)  x^1 


then  as  X  varies  between  0  and  1,  Xg^*^  is  continued  into  ,  Our  systc 

is  now 


I  <*<*>»  hU))  =  0  • 


As  before  we  use  (9)  and  Newton's  method,  with  X(X)  from 


*1  ai  .  af 
dx  =  ax  -  +  ax^T  h 


and  from  (16) 


*.  V  *,  *.  *, 

«_*  •  ’*«.*"•  '  •  '  •  ’’ «  * 
*  /  /  ’/  >  •  ■  •  •  • 

■  w.’  n-  m. .  »-  n-* 


,  *'«  •*,  »'.v 

•  •  4 


Pros  ($)  we  see  th«t  only  ^  depends  on  the  source  position.  Thus*  the 
Nx2  Matrix  df/dlg  has  all  eleaents  equal  to  aero  except  for  d#1/dx9  and 
d#t/dz#  in  the  first  row.  Since  our  program  uses  separate  subroutines  to 
wove  the  source  either  vertically.  (Fig.  3),  or  horisontally.  (Fig.  4).  one 
of  these  two  eleaents  is  zero  in  each  case.  The  subroutines  aay  be  combined 
to  aove  the  source  to  any  point,  though  this  method  is  particularly 
efficient  if  moving  on  a  grid.  Equation  (19)  rednces  to  the  vertical  or 
lateral  spacing  of  the  grid  points,  which  aay  differ,  another  advantage  of 
this  method. 


The  program  wet  written  in  FORTRAN  77  end  implemented  on  e  DEC  10  at 
the  Colorado  School  of  Nines.  There  are  three  entry  points  to  the  prograa. 
The  first  is  for  creating  ray  diagraas  like  those  in  this  paper.  This  part 
of  the  prograa  is  interactive  and  outputs  directly  to  the  terainal  or  to  the 
plotter.  The  user  supplies  a  subroutine  defining  the  interfaces,  their 
first  and  second  derivatives  and  the  interval  velocities.  Source  and 
receiver  positions  and  type  of  continuation  can  then  be  chosen. 

There  is  no  interaction  at  the  second  entry  point.  This  is  the  aass 
ray  tracing  procedure.  The  user  supplies  a  subroutine,  as  above,  and  also  a 
data  file  defining  the  region  of  interest i  the  location  of  this  region,  the 
nuaber  and  density  of  source  points  and  the  nuaber  of  receivers. 

Ve  start  by  deteraining  the  vertical  ray  in  a  horizontally  stratified 
aediua.  Ve  then  use  continuation  in  interfaces  to  obtain  the  correct  ray 
for  the  true  aodel.  The  interfaces  reaain  fixed  once  this  ray  has  been 
found.  The  region  of  interest  is  divided  into  a  grid  of  source  points  and 
we  aove  down  or  across  this  grid  using  source  continuation.  At  each  source 
position  we  use  continuation  in  receiver  location  to  find  the  rays  to  all 
receivers.  Each  raypath  is  stored  only  temporarily,  in  order  to  provide  a 
starting  point  in  the  solution  for  an  adjacent  ray.  After  each  ray  is 
found,  a  travel time  is  calculated  and  output,  so  that  for  each  source 
position  there  is  a  list  of  traveltimes  to  receivers.  Typically,  several 
hundred  thousand  traveltimes  will  be  found.  These  aay  be  used  later  in  this 
prograa  or  in  a  separate  inversion  scheme. 


In  each  of  our  continuation  schemes  no  start  by  setting  X  *  0,  AX  *  1, 
and  in  almost  all  eases  find  the  ray  solntion  in  one  step.  Typically,  two 
or  three  Newton  iterations  are  required  (and  never  more  than  four).  If  a 
solntion  is  not  found  within  four  iterations  we  return  to  the  known  raypath. 
set  AX  “  AX/2,  and  try  Newton's  method  again.  If  AX  falls  below  a  user 
defined  value  then  the  procedure  is  stopped.  Either  no  raypath  exists  or 
our  continuation  procedure  has  failed  to  converge.  It  is  sometimes  possible 
to  converge  on  a  solution  along  a  different  path  by  using  one  of  the  other 
continuation  procedures,  for  example,  by  returning  to  a  stratified  medium. 
Fawcett  discusses  such  cases  and  also  the  problem  of  bifurcation,  that  is, 
the  existence  of  many  solutions.  A  point  of  bifurcation  might  occur,  for 
example,  where  source  and  receiver  are  above  a  syncline.  Since  in  our 
scheme  the  source  is  at  depth,  we  avoid  many  such  situations. 

On  completion,  the  final  entry  point  will  output  a  migrated  depth 
section.  Each  point  at  depth  is  considered  a  possible  point  source.  The 
input  geophone  data  is  then  summed  over  curves  determined  by  the  traveltimes 
from  the  ray  tracing  procedure.  Only  at  points  along  an  interface  is  there 
significant  output  from  this  summation.  This  is  not  a  new  idea,  indeed  some 
of  the  earliest  computer  migration  schemes  used  existing  common  depth  point 
stacking  programs  that  summed  the  data  over  hyperbolae.  Our  medium  is  not 
of  constant  velocity  and  so  the  curves  are  no  longer  hyperbolae,  but  still 
the  traveltimes  serve  to  position  the  reflectors. 

It  is  a  "fast  and  dirty*  migration  and  we  make  no  claims  regarding  the 


amplitude  of  the  output 


(I)  Forward  Modeling 


Here  we  display  the  ▼ersatility  of  onr  technique.  Though  the  code  was 
designed  primarily  for  nse  in  inversion  it  can  be  applied  to  a  range  of 
forward  aodeling  probleas.  The  user  supplied  subroutine  defines  the  order 
of  interfaces  between  source  and  receiver  and  can  thus  specify  both 
reflections  and  refractions.  In  Figure  5  we  have  also  included  aultiples. 

By  placing  the  source  and  receiver  together  on  the  surface  we  force  a 
specular  reflection  froa  the  deepest  interface.  To  step  across  the  surface 
we  first  aove  the  receiver  and  then  the  source.  At  each  coincident 
source/ receiver  position  the  ray  is  plotted  or  the  traveltiae  is  calculated. 
(See  Fig.  6). 

In  a  siailar  fashion,  and  aotivated  by  our  interest  in  offset  probleas, 
we  first  separate  source  and  receiver  before  stepping  across  the  surface  to 
create  a  constant  offset  plot,  (Fig.  7). 

Vertical  seisaic  profiling  is  a  natural  application  for  vertical  source 
continuation,  (Fig.  8),  and  in  conjunction  with  receiver  continuation, 
(Fig.  9),  we  can  find  all  of  the  raypaths.  It  is  also  straightforward  to 
find  rays  reflected  froa  interfaces  below  the  receivers  in  the  well. 

Knowing  the  raypath  we  can  find  the  angle  at  which  it  strikes  an 
interface.  Also,  as  part  of  the  calculation  of  traveltiae,  we  found  the 


distance  travelled  within  each  layer.  That,  we  already  have  s one  factors  in 
the  calculation  of  aaplitnde  of  synthetic  data. 

(II)  Migration. 

The  only  synthetic  data  generator  available  was  for  beds  of  constant 
dip,  though  variations  in  thickness  and  interval  velocity  were  allowed.  We 
give  two  examples. 

Model  1. 

Figure  10  shows  three  interfaces  dipping  at  30  degrees.  The  interval 
velocities  are  5000,  6500,  10000  and  12000  feet/sec.  The  synthetic  data 
(Fig.  11)  is  for  fifty  receivers  spaced  at  100  feet.  The  frequency  of  the 
data  is  4  -  30  Hz.  As  we  night  expect  fron  the  "nigrator's  equation", 
(tan  (apparent  dip)  *  sin  (true  dip)),  the  slope  of  the  events  on  the  zero 
offset  tine  section  is  26.5  degrees.  Our  nigration  is  for  sone  region  of 
interest  above  which  the  positions  of  interfaces  and  the  interval  velocities 
are  known.  We  chose  the  region  bounded  by  dashed  lines  in  Figure  10.  The 
region  extends  fron  4100  feet  to  7050  feet  in  depth  and  is  4900  feet  across. 
The  horizontal  spacing  of  source  points  within  the  region  was  chosen  to  be 
100  feet  and  the  vertical  spacing  50  feet,  for  a  total  of  3000  source 
positions.  150,000  raypsths  were  found  in  tracing  to  all  fifty  receivers. 
C.P.U.  tine  to  detemine  these  raypsths  and  the  associated  traveltines  was 
15  ninutes  under  tinesharing  conditions.  The  sunnation  that  provided  the 
nigrated  output  shown  in  Figure  12  took  2  C.P.U.  ninutes.  The  interface 
has  been  positioned  accurately  and  its  dip  is  30  degrees.  More  receivers 
are  necessary  to  nap  the  interface  conpletely  within  the  region  chosen. 
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Model  2 


Here  we  have  increased  the  dip  to  75  degrees.  (Fig.  13).  The  synthetic 
data.  (Fig.  14),  was  generated  using  interval  velocities  10000,  12000,  12500 
and  14000  feet/sec,  with  100  receivers  spaced  at  150  feet.  The  frequency  of 
the  data  is  again  4-30  Hz.  The  region  of  interest,  shown  dashed  in 
Figure  13,  aeasures  900  feet  across  and  extends  froa  700  feet  to  4650  feet 
in  depth.  The  horizontal  source  spacing  was  chosen  to  be  100  feet  and  the 
vertical  spacing  50  feet.  The  output  of  the  aigration  is  shown  in  Figure 
15.  The  dip  is  74.5  degrees. 


Ve  have  extended  the  ideas  of  Perozzi  and  Fawcett  to  create  an 


algorithm  suitable  for  finding  large  nuabers  of  raypaths  and  traveltiaes. 
We  have  shown  that  the  technique  is  fast,  even  on  a  computer  of  liaited 
capacity.  Using  only  the  traveltiaes  froa  the  ray  tracing  procedure  we  have 
produced  highly  accurate  structural  inversions  for  two  specific  exaaples. 
Possible  applications  to  forward  aodeling  probleas  have  been  indicated. 
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We  describe  here  the  use  of  continuation  or  homotopy  procedures  in  an 
efficient  algorithm  for  the  calculation  of  raypaths  from  points  within  a 
region  at  depth  to  receivers  on  the  surface.  Our  work  follows  closely  that 
of  H.  B.  Keller  and  his  students  D.  J.  Perozzi  (1983)  and  J.  A.  Fawcett 
(1983).  Polynomials  define  the  interfaces  in  a  two-dimensional  piecewise 
constant  velocity  medium.  Starting  with  horizontal  layering  only,  the 
Interfaces  are  gradually  deformed  until  the  desired  earth  model  is  achieved. 
At  each  deformation  step  the  ray  equations  determined  by  Fermat's  Principle 
are  solved  using  Newton's  method  and  the  ray  from  the  first  source  position 
to  a  receiver  vertically  above  it  is  found. 

Next  we  employ  source  continuation,  moving  the  source  on  a  grid  within 
the  region  of  interest.  At  each  source  position  we  find  the  rays  to  all 
receivers  using  continuation  in  receiver  location. 

Knowing  the  raypaths,  it  is  straightforward  to  construct  a  table  of 
traveltimes.  These  traveltimes  alone  serve  to  position  reflectors  in  the 
subsurface. 


We  give  two  examples  of  migration  and  indicate  possible  applications  to 
forward  modeling. 
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